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We investigate the stability of the synchronization manifold in a ring and an open-ended chain 
(— 1 ■ of nearest neighbors coupled self- sustained systems, each self- sustained system consisting of multi- 

limit cycles van der Pol oscillators. Such model represents, for instance, coherent oscillations in 
biological systems through the case of an enzymatic-substrate reaction with ferroelectric behavior 
in brain waves model. The ring and open-ended chain of identical and non-identical oscillators are 
considered separately. By using the Master Stability Function approach (for the identical case) 
and the complex Kuramoto order parameter (for the non-identical case), we derive the stability 
boundaries of the synchronized manifold. We have found that synchronization occurs in a system 
of many coupled modified van der Pol oscillators and it is stable even in presence of a spread of 
parameters. 
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I. INTRODUCTION 



During the last decades, the emergence of collective dynamics in large networks of coupled units has been inves- 
tigated in disciplines such as physics, chemistry, biology and ecology. In particular, the effect of synchronization in 
systems of coupled oscillators nowadays provides a unifying framework for different phenomena observed in nature 

■ (for some reviews see [IHH]). Recently, complex networks have provided a challenging framework for the study of 
I synchronization of dynamical units, based on the interplay between complexity in the overall topology and local 

■ dynamical properties of the coupled units. A key problem is to assess conditions that guarantee the stability of the 
synchronous behavior for a network topology with some class of coupling configuration. To determine such conditions, 
many methods have been developed and among these, the so-called Master Stability Function (MSF) [J]. Indeed, 
the MSF approach was originally introduced for arrays of coupled oscillators and it has been later extended to 

H • the case of a complex networks of dynamical systems coupled with arbitrary topologies [5 8] . We will use the MSF 
JL" ' to investigate the stability boundaries of the synchronization manifold in a ring and an open-ended chain of self- 
. £^ \ sustained oscillators described by coupled multi-limit cycles van der Pol oscillators (MLC-vdPo) that might display 
birhythmicity, a mathematical model for instance suitable for some enzymatic substrate reaction with ferroelectric 
5— 1 ' behavior in a brain waves model @ , as well as in other biochemical models [l(| ■ 

Limit cycles are the fascinating objects in natural sciences where nonlinear kinetics due to feedback and co- 
operativity leads to instability of the fix point and evolves towards periodic sustained oscillations. Examples are 
abundant, with periods ranging from cardiac rhythms of seconds, glycolysis over minutes, circadian oscillations over 
the 24 hours, while epidemiological oscillations extend over the years |ld| - [l3j . Compared to the classical van der 
Pol system that possesses only one stable limit cycle, the MLC-vdPos has the advantage to better explain some 
biological processes, through its possibility of showing up multi limit cycles oscillatory states [lij. Thereby, we aim to 
understand how a large (but finite) number of interacting MLC-vdPos will behave collectively, given their individual 
dynamics and a corresponding coupling topology. We will investigate two oscillators coupling topology, a ring and an 
open-ended chain, both with nearest neighbors coupling. Indeed, such couplings are of paramount importance since 
they occur in many applications, both in natural and artificial systems, that involve the problem of coordination of 
multiple agents (circadian rhythm, contraction of coronary pacemaker cells, firing of memory neurons in the brain, 
superconducting Josephson junction arrays, design of oscillator circuits, sensor networks) [15l420l | . We will show that 
synchronization also occurs in the system of many coupled modified van der Pol oscillators and it is stable even in 
presence of a spread of parameters. 

The organization of the paper is the following: Section II deals with the description of the model under consideration. 
In Section III, the stability boundaries of the synchronized states in a ring of identical MLC-vdPos are retrieved through 
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the MSF approach, and numerical simulations are used to validate and complement the results. Section IV addresses 
the ring of non-identical MLC-vdPos, while Section V considers the stability of synchronization manifold in a chain 
of open-ended nearest neighbors coupled identical and non-identical self-sustained oscillators. Finally, Section VI is 
devoted to the conclusions. 



II. THE SELF-SUSTAINED MODIFIED VAN DER POL OSCILLATOR 

The modified van der Pol oscillator is described by the following nonlinear differential equation (in non-dimensional 
form) 

x c — /Lt(l — x 2 c + ax* — fix®)x c + x c = 0. (1) 

(Overdots stand for the derivative with respect to time). Such model, that exhibits an extremely rich bifurcation 
behavior, was proposed by Kaiser [l4j |. and describes many dynamical systems, ranging from physics, biochemistry to 
engineering. When employed to model biochemical systems, namely the enzymatic-substrate reaction, x c represents 
the rate of change of the number of excited enzyme molecules and x c in Eq. (1) is proportional to the population of 
enzyme molecules in the excited polar state. The quantities a and j3 are positive parameters which measure the degree 
of tendency of the system to a ferroelectric instability compared to its electric resistance, while fi is the parameter 
that tunes nonlineari ty [2 l|. The nonlinear dynamics and the synchronization process of two such systems have been 
investigated recently |2ll.l22j] , while Ref . [23| considers its dynamics and active control to find that chaos can be tamed 
for an appropriate choice of the coupling parameters. Coupled van der Pol oscillators have been considered in (23 - [26| . 
Depending to the values of the parameters (3 and a, Eq. (jTJ) can lead to one or three limit cycles. When three limit 
cycles are obtained, two of them are stable and one is unstable, a condition for birhythmicity [lOj: The unstable limit 
cycle represents the separatrix between the basins of attraction of the two stable limit cycles (see |ll|). We will avoid 
such region for simplicity, so the uncoupled elements are monorhythmical. We show in Fig.l the region of existence 
of birhythmicity in the two parameter phase space (f3-a) [2l|, |22j . 

III. THE RING OF IDENTICAL SELF-SUSTAINED OSCILLATORS 

A. The model 

When the coupling is realized through a ring topology, i.e. with periodic boundary conditions, the model is described 
by the following set of dimensionless nonlinear differential equations: 

'i\ — fi(l — x\ + ax l — (3x\)ii + X\ = K{x% — 1x\ + xn), 

x v - /x(l - xl + ax v - (3xl)x u + x v = K(x u+ i ~ 2x v + x„_i), v = 2,3, N - 1, 

xn — Ml - 4 + ax N ~ P x n) x n + xn — K(xi — 2x N + (2) 

Here, K stands for the diffusive (nearest-neighbor) coupling strength. The physical meaning of the variables depends 
on the nature of the system described by such a ring. For instance, let us consider the case of a ring of immobilized self- 
sustained enzymes that minimize the cost of production by permitting repeated use of the enzymes and substantially 
increases the stability of the enzyme reactions themselves [27j. We propose to describe the coupling between cells 
by means of the diffusion of the solute concentration, so K is proportional to the difference in the concentration. 
Since within each self-sustained system there are promoters (positive direction of solute flow) and inhibitors (negative 
direction of solute flow) that determine the direction of the solute flow, it follows that the coupling coefficient K can 
be either negative or positive. 

The system of diffusive coupling occurs in many complex systems. For instance, the ring of multi-limit cycles 
oscillators is of interest to construct hypotheses on the in vivo enzymes behavior when they possess more than one 
stable limit cycle. A ring of the MLC-vdPos can therefore serve as a simple model of a biological oscillator. The 
set of equations © can also be exploited in electronics engineering as a network of parallel microwaves oscillators 
28l [29j : Multiple limit cycles can be obtained from a classical van der Pol oscillator with hysteresis in the inductance 
30]. Since we are interested in the synchronization manifold, the following subsection deals with the stability of the 
synchronous states in the ring of mutually coupled identical self-sustained oscillators. 
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B. Stability of the synchronization manifold 

We aim now to determine the stability of the synchronous state in the system (2) of self-sustained oscillator requiring 
that each of the perturbed trajectories returns to its original limit cycle. We are therefore interested in the bifurcations 
from the state that resides on a synchronization manifold denoted by M. ={(xi,yi) = (£2,2/2) = ••■ = (xniVn)}- In 
order to rewrite Eqs. ([2]) as set of flows, we have introduced a new variable y v = x v \ with such new variables the 
corresponding dynamic equations read: 



Vv = mC 1 ~ x l + ax t ~ P%l)yv - Xy + K(x u+1 - 2x v + x v -i), 



1,2,...,N. 



(3) 



To verify the stability of the synchronization manifold M. we make use of the MSF approach 0, Thereby, let 
X 1 be the two-dimensional vector of the dynamical variables of the i th unit, H : R 2 — > R 2 an arbitrary function 
describing the coupling between each unit variables. Thus, the dynamics of the i th unit is rewritten as a function of 
the 2 x N column vector state X 1 as 



N 



X i = F(X ; ) + KJ2 GyH(Xl), 



1,2,..., N, 



(4) 



where X.*=[ Xl , Vi } T , F(X ; )= [ yi ,fx{l 



ax; 



E = 



" Xi\ 



1 



, and the function H is defined through the matrix 



by H(X') = EX 1 . Gij £ R are the elements of the N x N symmetry connectivity matrix G defined as 



G 



/ -2 1 
1-2 1 
1 -2 



V 1 



1 



1 \ 




-2 J 



Since the ring is made of identical self-sustained oscillators, the evolution function -F(Xi) in Eq. (4) is the same for 
all ring node. This ensures the existence of an invariant set (xjft), yi(t)) = (x s (t), y s (t)), Vi, representing the complete 
synchronization manifold Ai . Following the MSF scheme [l|, Q, E| , 
be determined by letting 



the stability of the resulting dynamical states can 



x v (t) = 8x v {t) + x s (t), 
y v {t)=8y u {t)+y s {t), 

and linearizing equations (j4]) around the periodic limit-cycle state (x s ,y s ). This approach leads to the following 
equation: 



5±(t) = [1 N (g) JF(X S ) + i^G (g) 3U{X,)]6X(t), (5) 

where (^) stands for the direct product between matrices. J denotes the Jacobian operator and the 2 x N column 
vector SK. l (t)= (8xi(t), 8yi(t)) is the deviation of the i th vector state from the synchronization manifold. We have 
used the definitions H(X i ) = EX', and JH = E. 

A necessary condition for stability of the synchronization manifold is that the set of (N — 1) x 2 Lyapunov exponents 
that corresponds to phase space directions transverse to the 2-dimensional hyperplane X*(t) = X s (f) should be 
entirely made of negative values [3, [3l|. The arbitrary state JX(i) can be written as <5X(£) = J^iLi v i with 
= te,i(i))^i,2(i))- If one applies [vj] T (7, and [vj\ are the set of real eigenvalues and the associated orthonormal 
eigenvectors of the matrix G) to the left side of each equation in ([5]), one finally obtains the following set of N 
variational equations 



- [JF(X S ) + X 7fe JH(X s )]£ fc (i), k = 0, 1, 2, N —1. 



(6) 
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The eigenvalues 7^ of G are given by 7^ = — 4sin 2 (7rfc/iV) 0, [Hj]. Let us remark that the mode k = is the "slower" 
or uniform mode, while k = N — 1 is the "faster" or more rapidly oscillating (in the momentum space) mode. Each 
mode k in Eq. ([6|) corresponds to a set of 2 conditional Lyapunov exponent Xj, (j = 1, 2) along the eigenmode related 
to a specific eigenvalue 7^. Eqs. ([5]) enable to compute the maximum Lyapunov exponent \™ ax of each mode k as 
a function of the coupling parameter K. In the following, we will use the parametrical behavior of the largest of 
such exponents A(if), called MSF, to determine the stability boundaries of the synchronization states in the ring of 
coupled self-sustained systems. 

C. Numerical simulations 

We now address the numerical computation of the Lyapunov exponents and then the MSF A, by solving the 
equations of motion (1) and the variational equations (6), with a fourth-order Runge-Kutta algorithm. The parameters 
used are representative of the monorhythmical solutions: /1 = 0.1, and in (i): a — f3 — 0.1, while in case (ii) we set 
a = 0.114, /3 = 0.005. 

When the coupling coefficient is turned off (i.e. K = 0), the self-sustained systems are uncoupled and the corre- 
sponding MSF is A(K = 0) = 0. For K 7^ 0, and for a given ring size N, the MSF A(K) will enable us to derive the 
range of the coupling parameter K in which the transverse Fourier modes are stable, and therefore each oscillator 
of the ring exhibits the same dynamical states. The behavior of the MSF as a function of the coupling K and the 
number of oscillators N is plotted in Fig. 2. Starting from large negative K values, the first domain is the unstable 
synchronization domain (Djjs) in which there is some positive transverse Lyapunov exponent. The boundaries slowly 
change with the number N of oscillators (see Table 1); when synchronization is not observed, all the modes are on 
the transverse manifold where variations transverse to the synchronization manifold do not decay; all or some of the 
transverse Lyapunov exponents are positive, i.e. \ max > 0. 
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Domains of coefficient K for unstable synchronization (Djjs) 
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Table 1: Range of coupling values for unstable synchronous state in the ring of oscillators with one limit-cycle. Parameters of the system are 

n = 0.1, a = 0.1, p = 0.1. 

In other words, the Fourier modes increase continuously or possess a bounded oscillatory behavior and one never 
observes stable synchronization in the ring for such a choice of the coupling strength. It should be underlined 
that the region Djjs where unstable synchronization is found, can be divided in two sub-domains: the first sub- 
domain corresponds to the non-synchronization phenomena, while the second sub-domain corresponds to partial 
synchronization. For example for N = 4, in the first sub-domain (K g] — oo; — 0.25[), any perturbed trajectory leads 
the oscillators to continuously drift away from their original limit cycles. In the second sub-domain, (K e] — 0.02; 0[) 
although perfect synchronization is not observed, there is some tendency towards co-operative behavior, as will be 
shown in Sect. IV by means of the Kuramoto order parameter. 

It appears that as the coupling coefficient K increases from — oo, where the ring of coupled self-sustained systems 
is on the unstable synchronized states (and therefore all Lyapunov exponents are positive), the mode k = 1 is the 
first one to move from the unstable to the stable domain. On the other hand, k = N/2 is the last mode that leaves 
from the unstable domain to enter the stable state. At this point, the ring is on the stable synchronized state; when 
the coupling coefficient K is further increased, the mode k = 1 is the first mode to move from the stable domain 
to the unstable one, and the mode k = N/2 is the last one. For the case of the positive coupling coefficient K, the 
reverse situation is observed: the k = N/2 mode is the first to become stable. Also in the case (ii) with a = 0.114 
and f3 = 0.005, is found a similar behavior. 

Since from a biological point of view, negative and positive values of the coupling strength refer to inhibitory and 
excitatory process, respectively, one can conclude the following: The instability of a synchronization process within a 
system of inhibitory (K < 0) cells is driven by the slowest, or lower k, mode. On the other hand, when a system is 
made of excitatory (K > 0) cells, the instability of the process is initiated through the fastest, or higher k, mode. 

Fig. 3 shows the stability diagram in the (N, K) plane, i.e. it displays the border between the two main dynamical 
states. We have found with the MSF method that the number of units strongly affects the stability boundaries of 
the synchronization process in the ring. In particular, when N > 30, there is no domain of stable synchronization for 
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negative values of the coupling coefficient K. ft is important to notice that Eqs. ([6]) have been linearized around the 
stable orbit. In order to check the validity of the above semi-analytical results, we have numerically solved Eqs. ©. 
Numerically, the stable synchronization state is achieved if all the oscillators of the system are in synchrony. This 
means that the synchronized state defined by the N constraints (x\,y\) = (0:2,1/2) — (^3>2/3) = ■■■ = ( x n,Un) (within 
some precision or tolerance h) is in the synchronization manifold and therefore the stability of the synchronization 
process is ensured. The numerical synchronization criterion is defined as the sum of the absolute deviations from each 
oscillator respect to all other, and dividing such sum for the number of terms: 



It should be underlined that this numerical synchronization criterion gives an insight about perfect synchronization 
in the ring but does not tell anything about other types of dynamical behavior related to stable synchronization. 
In some regions where stable synchronization is expected from MSF analysis, one can instead observe unstable 
synchronization in the full model (2). The discrepancy depends upon the choice of initial conditions and the number 
of the self-sustained oscillators which are taken into account during the process. The neglected nonlinear terms 
give rise to nonlinear contributions that apparently increase with increasing the coupling K, thus explaining the 
discrepancy between the MSF method and the numerical simulations. From the criterion ([7]), we have obtained some 
results that depend on the number N. In particular we have noticed that increasing the number of oscillators the 
domain of stability of all modes shrinks. For instance, increasing the number of oscillators from N = 6 to N = 8 in 
Fig. 3 (ii) the Lyapunov exponents associated to each mode become positive for a lower value of the coupling constant 
K: the first mode becomes unstable for K > -0.014 when N = 6, for K > -0.019 for TV = 7, and for K > -0.024 
for N — 8. An analogous behavior is observed for the other modes. 

Both full synchronized states and no synchronization can be found in Fig. 3, where we also display the comparison 
between analytical and numerical results. One also finds some domains obtained only from numerical simulation that 
do not match with the MSF stability condition. It is important to notice that this gap depends upon the choice of the 
precision or tolerance h. In fact, given the accuracy of the numerical integration, if h is too small (h << 10~ 6 ) one 
would never observe synchronization, while for too a large h (h >> 10 -6 ) synchronization would appear as a spurious 
effect, so it is necessary to determine an appropriate value of the parameter h as a function of the accuracy of the 
numerical scheme, a value that we have estimated for our scheme to be h ~ 10 -6 . To circumvent such difficulty (and 
arbitrariness) we will suggest in Sect. IV another tool to detect synchronization, the Kuramoto order parameter. 

Figs. 4 and 5 show space-time-amplitude diagrams that display the behavior for two values of the coupling parameter 
K in the unstable or no-synchronization (US) and stable synchronization (SS) domains. For the unstable domain, 
Fig. 4 shows the time evolution for (N,K) = (30,-0.2) and (N,K) = (50,-0.05), where no synchronization is 
found in the ring. Considering the domain of stable synchronization, we display in Fig. 5 for (N, K) — (20, 0.6) and 
(N,K) = (50, 1) the process of entrainment of the oscillators in the ring. 



The understanding of the collective dynamics of synchronization processes occurring in living organisms is a 
paramount task of formidable difficulty. One reason is that in a realistic description of such systems, even within the 
same species, oscillators are made of non identical elements. Although not explicitly considered in this work, real 
systems are also always subject to noise in the form of fluctuations associated with dissipation, as well as in the form 
of random force due to the external environment. We will, however, focus our attention on quenched disorder, or 
non identical oscillators that can be modelled with a spread of the characteristic parameters. For example, a ring 
of nonidentical oscillators is useful to investigate calcium dynamics in pancreatic acinar cells (33| as well as other 
complex networks [34]. In the case where Eq.@ is used to model biological systems of immobilized enzymes, the 
interpretation of the parameters /1, a and /3 requires that they can vary from site to site, according to the tendency 
degree of the oscillator to possess a ferroelectric behavior and also to the conductivity of the medium [2l[ . This leads 
to the possibility to have nonidentical biological oscillators, or to a ring of nonidentical coupled systems described by 
the following equations, where a and /3 are constant, while \i is not: 




(7) 



IV. THE RING OF NONIDENTICAL COUPLED OSCILLATORS 



x\ — /ii(l — x\ + ax\ — f3x\)xi + w\x\ = K{x% — 2x\ + xn), 

x v - jUr/(l - xl + ctx\, - pxl)x u + wlxy = K(x v+ i - 1x v + x„_i), 

'xn - MJv(l - x n + ax tr - Px%)±n + w 2 n xn = K(xi - 2x N + zjv-i)- 



v 



N-l, 



(8) 
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We recall that \i v are positive coefficients; they measure the dissipative strength and are smaller than unity. w v is 
the frequency of the associate linear oscillator (/x„ = 0) and depends upon the site. The nonlinear terms play the 
role of amplitude-dependent dissipation and provide a self-sustaining mechanism for the perpetual oscillation. For 
the original biological system [2l|, |35| , fi u can be expressed as function of the parameters of the systems as follows 



K — u 18 2 

= = —aWovK. , (9) 

w ov 5 

where w ov are the frequencies of the periodic enzyme-substrate reaction and can be determined by the recombination 
and attraction coefficients ^Yvi^v as w ov — VtvCT (21j - £ v is the decay rate of each excited enzymes to the ground (or 
weakly polar) state and 7^ the range attraction of the substrate particles due to the autocatalytic reactions. From 
Frolich ideas, we may suppose that in large regions of the system of proteins, substrates, ions and structured water 
are activated by the chemical energy available from substrate enzyme reactions [36| . Thus, chemical oscillations in the 
number of substrate and activated enzyme molecules with a very low frequency w ov might be carried out around the 
equilibrium state [35} ■ k (obtained through a nonlinear dielectric contribution) is the coefficient proportional to the 
macroscopic polarization and the time dependent number of the excited enzyme molecules [35| . It can be chosen in the 
interval 1 < n < 5. a is viewed as a coefficient of relaxation term of electric resistances against the systems tendency 
to become ferroelectric. We assume that the frequency w ou is chosen as: w ov = 1 + Au; (C„ — i) (( u is a uniform 
random variable less than the unity and Aw a disorder parameter). We also assume that the natural frequencies 
w v are uniformly distributed around 1: w v = 1 + Au>(£„ — ^). For simplicity, to avoid a new parameter, we take 
the spread of the recombination frequencies Aw and of the natural frequencies Aw to be the same: Aw a = Aw. 
According to the analytical expression the amplitudes A and frequencies fl u of the limit-cycles established in Ref . [2l| . 
the periodic solution for each oscillator is approximately described by 

x v (t) = Acosfl,yt (10) 

where the amplitude A and the frequencies fi„ for the MLC-vdPo (P) are derived through the following equations: 

f^-^ + I^-l = 0, 11(a) 
, fj, 2 f 1580/3 „ 12 738a/3 10 , 72a 2 + 309/? , 8 , 64a - 219/3 , 6 

\l v -W v + — <— — — A - nnno . A +[ — )A -( — — )A 



w 



v \ 393216 99024 v 768 ' y 6144 



16a + 3 4 _£ 2 1 3 
v 384 ; 64 J u ; v 1 

We note that only the frequencies Vl u are modified by the disorder parameter Aw , while the amplitude A of the 
orbit does not change. Thus, each self-sustained system in the ring exhibits a different frequency, the amplitude, 
see Eq. ([TIT a)), only depends upon the a and /3 parameters, which are fixed here. The self-sustained oscillators 
with nonidentical /x„ consist therefore of oscillators with nonidentical natural frequencies £l v , that are a function 
of the physical parameters /x„, a, /3, while the amplitude A of the limit-cycles is unchanged [21]. Therefore in the 
uncoupled limit approximated by the harmonic oscillations, see Eq.(10), each oscillator is completely described by 
a phase S v = cos -1 (x v (t)/A) and amounts to a rotator. It is therefore tempting to monitor the rotators through 
the phases 5„ analogously to the Kuramoto model l37l-[39l . a paradigm for coupled rotators and already employed 
to investigate globally coupled van der Pol oscillators [40|. Synchronization of the nonidentical oscillators (8) can be 
revealed by means of the Kuramoto order parameter R defined as (37T439T ] : 



i?=lf>^-. (12) 



For a completely disordered steady-state one gets (for very large N) R ~ 0, while R ~ 1 corresponds to phase synchro- 
nization. The advantage to employ such parameter R is that it is possible to measure the degree of synchronization, 
or, roughly speaking, the fraction of synchronized oscillators in the intermediate cases of partial synchronization, for 
which the parameter R will reach a value < R < 1. So the Kuramoto order parameter R is capable to describe also 
states where condition ([7]) is not satisfied, but nevertheless the oscillators exhibit some tendency to behave coherently. 
The interest to determine the tendency of nonidentical systems to behave synchronously is almost ubiquitous, from 



laser arrays [41| to people walking on a bridge [42j, or Josephson arrays [43J. For instance in the case of Josephson 



oscillators it has been speculated that the analysis of the Kuramoto order parameter (|T2|) is useful to predict the 
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fraction of coherently working electronic elements, and therefore the efficiency of the devices [lg, |44|, |45[ . Here we will 
use the Kuramoto order parameter in the same spirit: as a diagnostic tool to measure the capability of the oscillators 
to synchronize in spite of the differences in the single oscillators frequencies. In order to look at the qualitative picture 
of the dynamical states of oscillators under the influence of diffusive coupling, we have plotted in Fig. 6 the Kuramoto 
order parameter as a function of the coupling parameter K for several different finite number N of oscillators, as 
analyzed for instance in Ref. [46[. From Fig. 6, it can be noticed that full synchronization occurs mainly for K > 0, as 
expected from MSF analysis, although the Kuramoto order parameter allows an analysis of the partial synchronization 
( corresponding to a state - analogous to clustering - where the maximum Lyapunov exponent is positive, but not all 
modes are unstable 47]) that occurs in the region where full synchronization is unstable, see Sect. III. Therefore, this 
finding suggests that in the case of non-identical oscillators, excitatory coupling is more suitable to achieve coherent 
oscillations when compared to the inhibitory coupling (K < 0). It can be notice that synchronization in presence of 
disorder is not achieved for negative K values as we show on Fig. 7 for an enlargement around low K values. From the 
Kuramoto order parameter (TT2|) . the following results have been obtained for N = 6, N = 10, N — 20 and N = 30. 
When the ring is composed of N = 6 nonidentical coupled oscillators, the phase synchronization state is achieved for 
K > 1.6. As the number of nonidentical oscillators in the ring increases, one still observes two main dynamical states 
but for different ranges of K, and in general the degree of synchronization decreases when the number of oscillators is 
increased. In the case of a ring is composed of N = 10 and 20 nonidentical oscillators, full synchronization is obtained 
approximately in the range K > 3.6 and K > 10.5, respectively. For N = 30 oscillators in the ring, the oscillators are 
fully synchronized only for K > 26.0. 

It is also interesting to compare the result of the Kuramoto order parameter analysis in Fig. 7 with the MSF analysis 
of Fig. 3. In this enlargement it is clear that disorder (Aw — 0.05) induces an enlargement of the desynchronization 
region around K = 0. For instance for N — 10 while the uniform oscillators are not synchronized only in the 
narrow region —0.028 < K < 0.005, after which there is a range of perfect synchronization (see also Table 1), from 
the simulation of disordered oscillators in Fig. 7 it is evident that there is a much more extended region of partial 
synchronization. Nevertheless it is clear that the qualitative behavior is the same in Fig. 3 and Fig. 6: the region 
unfavorable to synchronization around low K values clearly increases increasing the number N of oscillators in both 
the ordered and disordered cases. 

The behaviors of the space-time-amplitude in the partially synchronized state and in the fully synchronous motion 
are similar to the patterns shown in Figs. 4 and 5. The ring of nonidentical oscillators capability to give rise to 
synchronization is represented in Table 2. 
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Domains of K for synchronization (ii) 


2 


[0.4;+oo [ 


[0.2;+oo [ 


3 


[0.6;+oo [ 


[0.3;+oo [ 


4 


[0.7;+oo [ 


[0.6;+oo [ 


5 


[1.0;+oo [ 


[0.7;+oo [ 


6 


[1.3;+oo [ 


[0.8;+oo [ 


8 


[2.4;+oo [ 


[1.3;+oo [ 


10 


[3.6;+oo [ 


[1.6;+oo [ 


15 


[7.2;+oo [ 


[4.1;+oo [ 


20 


[12.8;+oc [ 


[7.6;+oo [ 


25 


[19.4;+oc [ 


[ 10.6: + oo [ 


30 


[26.0;+oc [ 


[14.8;+oc [ 



Table 2: Domains of full (R > 0.998J synchronization in the ring of nearest neighbors coupled nonidentical oscillators. Parameters of the 
system are (i) — 0.1, a — 0.1, /3 — 0.1; (ii) — 0.1, a — 0.114, /3 — 0.005. The spread of the frequencies is Aw — 0.05. 

Figure 8 shows the stability diagram in the (N, K) plane and reveals the boundary between phase synchronization 
(defined as R > 0.998) and disordered steady-state. These curves also show that the dynamical states depend on 
the ring size. Analyzing the effects of the disorder parameter Aw (that controls the spread of the frequencies of 
the uncoupled oscillators) on various dynamical states in the ring, one finds that such disorder parameter can also 
induce de-synchronization, i.e. even preparing the oscillators with identical initial conditions, the spread induces a 
disordered state. 
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V. AN OPEN-ENDED CHAIN OF NEAREST NEIGHBORS COUPLED OSCILLATORS 
A. A chain of open-ended nearest neighbors coupled identical oscillators 



N 


Domains of K for stable synchronization (i) 


Domains of K for stable synchronization (ii) 


2 


]-0.5;-0.01[U]0.014; + oo [ 


]-0.5;-0.01[U]0.009;+oo [ 


3 


]-0.34;-0.02[U]0.02; + oo [ 


]-0.334;-0.019[U]0.019; + oo [ 


4 


]-0.3;-0.03[U]0.04;+oo [ 


]-0.3;-0.03[U]0.03; + oo [ 


5 


]-0.28;-0.05[U]0.07;+oo [ 


]-0.28;-0.05[U]0.05;+oc [ 


6 


]-0.27;-0.07[u]0.1;+oo [ 


]-0.27;-0.07[U]0.07;+oc [ 


8 


]-0.26;-0.13[U]0.18:+oo [ 


]-0.26;-0.13[U]0.12;+oo [ 


10 


]-0.26;-0.21[u]0.29:+oo [ 


]-0.26;-0.21[u]0.19;+oo [ 


20 


]1.16;+oo [ 


]0.77: + oo [ 


30 


]2.61; + oc [ 


]1.74; + oo [ 


40 


]4.66;+oc [ 


] 3.1: + oo [ 


50 


]7.28;+co [ 


]4.85:+oo [ 



Table 3: Domains of stable synchronization in the open-ended chain of nearest neighbors coupled identical oscillators obtained with the MSF 
approach. Parameters of the system are (i): — 0.1, a — 0.1, (3 — 0.1, (ii) p. — 0.1, a — 0.114,/9 — 0.005. 

Not closed chains of coupled nonlinear systems, corresponding to open-ended boundary condition, are involved in 
many natural processes such as the swimming motion of organisms [48j and waves synchronization that occur during 
sensory processing in the cortex [49L [50|. In the case of small intestinal muscle, it has proved useful to simulate and 
attempt an analysis of long chains of oscillators [H| ■ There is in fact physiological evidence that a complete organ such 
as small intestine, comprises a very large number of smooth muscle cells organized to form self-oscillatory segments. 
Also, the electrical waveforms found in the canine gastrointestinal tract are very non-sinusoidal in nature, whereas 
research work realized in the De par tment of Surgery, University of Sheffield has indicated that the human duodenal 
signals are nearly sinusoidal [52-54], Consequently, the nonlinear parameter fi can be assumed small to model the 
duodenal oscillator cells. In fact in the autonomous regime, the final limit cycle state of the van der Pol oscillator is 
nearest a sinusoidal behavior for small values of fi. In contrast, for large [i it develops relaxation oscillations [55l. |56|. 
In the open-ended chain coupling configuration, one assumes that the 1 th and N th oscillators are not interconnected. 
In this case, the differential equations of motion are defined as 



x\ — — x\ + ax\ — /3xf)ii + Xi = K(x2 — xi), 
Xi — — xf + ax\ — fix\)xi + Xi — K(xi + i — 2xi + x 
x N — fi(l — x 2 N + ax% — flx%)x N + x N — K(xn- 



i = 2,3, ...,N- 1 



x N ). 



(13) 



Our purpose is to recover the essential features of the stability investigation (employing the same methods of the 
section III and IV) to identify some general properties of an open-ended chain of nearest neighbors coupled self- 
sustained systems. The stability of the synchronization process is again found following the MSF approach. In the 
present case, the coupling matrix is 



G 



OE 



/-I 
1 





\ 





V 



-1/ 



N 


Domains of K for stable synchronization in the ring 


Domains of K for stable synchronization in the open-ended chain 


10 


[-0.24;-0.03]U[0.06;+oo [ 


]-0.25;-0.22[u]0.3:+oo [ 


20 


[-0.24;-0.12]u[0.22;+co [ 


]1.17;+oo [ 


30 


[0.49;+oo [ 


[2.62;+oo [ 


40 


[0.87;+oc [ 


[ 4.66; + oo [ 


50 


[1.35;+oc [ 


[7.28;+oo [ 



Table 4 : Comparison between stability domains in the ring and the open-ended chain of identical oscillators obtained with the MSF approach, 

case (i): fj, = 0.1, a = 0.1, (3 = 0.1. 
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The coupling matrix Goe also obeys the zero sum condition, and the synchronization manifold A4 is an invariant 
set. Therefore, stability of the synchronous state reduces to find the behaviors of the system dynamical properties 
along directions in phase that are transverse to the synchronization manifold. The stability of the resulting dynamical 
states can be determined considering the set of variational equations (6). In Rcf. [32], it was showed that the matrix 
Goe can be diagonalized in a manner similar to the shift-invariant case using a discrete Fourier transform with N 
replaced by 27V, and the eigenvalues are jk — —4 sin 2 (irk/2N) , for k = 0, 1, N — 1. This is similar to the case 
of diffusively coupled self-sustained systems, except that there are no degenerate modes and the highest wavelength 
corresponds to k = N — 1. Nonetheless, because of the dependence of the eigenvalues on TV, we will see the same 
dynamical states as before. The stability boundaries of the synchronization derived through the MSF are shown in 
the Table 3: 

The dependence of the size of a chain also apparent in Fig. 9, where the MSF is plotted in the plan (K, TV) for the 
open-ended chain. It appears that, analogously to the results reported in Fig. 2, as the size of a chain increases, the 
domain of stable synchronization reduces more quickly than in the case of the closed end model, see Table 4. 

In Fig. 10 we plot the stability diagram in the (TV, K) plane for the open-ended chain of coupled self-sustained 
oscillators, obtained through the MSF and a direct numerical simulations of the differential equations (|13[) . Comparison 
with the case of the ring, (see Fig. 3), reveals that the agreement between the results obtained from the MSF and 
those of a direct numerical simulation is very good. It also indicates that the region of stable synchronization is 
smaller then in the case of the ring, as reported in Tables 4 and 5. 



N 


Domains of K for stable synchronization in the ring 


Domains of K for stable synchronization in the open-ended chain 


10 


[-0.24;-0.04]U[0.04;+oo [ 


]-0.25;-0.21[u]0.2;+oo [ 


20 


[-0.24;-0.12]u[0.13;+co [ 


]0.78;+oo [ 


30 


[0.29;+oo [ 


[1.75;+oo [ 


40 


[0.51;+oo [ 


[ 3.11: + oo [ 


50 


[0.79;+oc [ 


[4.86;+oo [ 



Table 5: Comparison between stability domains in the ring and the open-ended chain of coupled identical oscillators obtained with the MSF 

approach, case (ii): n — 0.1, a — 0.114,/3 — 0.005.. 

Finally, let us add that the space profiles and the time dependent amplitudes in both cases of stable and unstable 
(or partially synchronized) states are very similar to those shown in Figs. 4 and 5, and therefore are not shown. 

B. A chain of open-ended nearest neighbors coupled non-identical oscillators 

In the case of an open-ended chain of nearest neighbors coupled non-identical oscillators, the chain's dynamics is 
given by the following set of equations: 

'i\ — /ii(l — x\ + ax\ — flx\)x\ + w\x\ = K{x2 — xi), 

ii - fii(l — x 2 + ax\ - fix1)±i + wfxi = K(x i+ i - 2xi + Xi-i), i = 2, 3, TV - 1, 

xn - Mjv(1 ~ x li + ax A N - f3x%)x N + w 2 N x N = K(x N ^i - x N ). (14) 

Our choice of considering non-identical coefficients is justified by the parameter spread encountered in real systems, 
as discussed in section IV for the ring of nearest neighbors coupled non identical self-sustained oscillators. The 
set of equations (TT4"| is numerically integrated to compute the Kuramoto order parameter R, Eq. (|12p . versus the 
coupling coefficient K, to quantify the fraction of synchronized oscillators dynamics in the chain of open-ended nearest 
neighbors coupled non-identical oscillators. 

Fig. 11 shows the variation of the Kuramoto order parameter R versus the coupling coefficient K for TV = 6, 
TV = 10, TV = 20 and TV = 30. Phase synchronization is found in the open-ended chain and the stability boundary 
depends on N. However, as the number of oscillators TV increases, the open ended chain is very sensitive to the initial 
conditions and this dependence more pronounced when the coupling strength K becomes larger. Fig. 12 shows the 
stability diagram for an open-ended chain of coupled nonidentical oscillators in the (K, TV) plane, and it is qualitatively 
very similar to Fig. 8. The chain of coupled nonidentical oscillators capability to give rise to synchronization (defined 
again as R > 0.998) is appears in Fig. 12 represented in Table 6 for both the sets of parameters. 

VI. CONCLUSIONS 

We have studied some criteria under which the synchronization manifold is stable for a ring and an open-ended 
chain of nearest neighbors coupled van der Pol-like oscillators. By using the Master Stability Function, the stability 
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boundaries of the main synchronized states have been investigated and the obtained results have been complemented 
by numerical simulations. Thereby, the following findings have been captured: the threshold of the coupling strength 
from which the system displays a stable synchronized behavior is different for the cases of two sets of parameters, 
as depicted in Fig. 2. Moreover, as far as the diffusive coupling is concerned, we have found that with just the first 
and last mode's behavior one can characterize the stability of the system (as far as full synchronization is considered, 
i.e. neglecting clustered states). The stability boundaries of the synchronization have been derived as a function of 
the number of oscillators and of the coupling strength. We have found that increasing the number of oscillators, the 
synchronized states are unfavored, while, quite obviously, a larger coupling constant increases their stability. The non 
obvious result is the presence of a stable synchronization manifold even for negative values of the coupling K, i.e. 
for repulsive interaction. We have also found that, when the size of the ring and the coupling coefficient increases, 
the domains of instability become very large and the gap between analytical and numerical analysis becomes more 
relevant. 

In the case of non identical oscillators, we have performed numerical simulations to compute the fraction of entrained 
oscillators by means of the Kuramoto order parameter. Varying the coupling strength, we have retrieved the region 
of unstable and stable synchronization states of the ring and an open ended chain of non-identical nearest neighbors 
coupled self-sustained oscillators. In particular, we have found that synchronization is robust enough to survive also 
in presence of some amount of disorder. However, the process depends upon the number of oscillators and the amount 
of disorder. 

Research might be extended in several directions. For example we think that an extension of the analytic treatment 
to find stability of synchronous states in the ring of identical and non identical coupled oscillators in the presence 
of noise is an interesting task which can be tackled also using the Kuramoto order parameter. Also, we have re- 
stricted our research to a region of the parameters where only monorhythmical states are present. An interesting 
question is therefore the influence of two stable orbits on the synchronization properties in presence of disorder. 
Moreover, the spontaneous transition frequencies in the ring with external random excitation, such as noise, can be 
also experimentally observed if one can resolve in time enzymatic reactions or electronic circuits. 
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FIG. 2: Variation of the Master Stability Function (MSF) A versus K and N for a ring of diffusive coupling oscillators, (i) 
H = P = a = 0.1 and (n) fj, = 0.1; f3 = 0.005, a = 0.114. 
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FIG. 3: Stability diagram in the (N, K) plane of a ring of diffusivly coupled oscillators. The shaded area denotes the region of 
stable synchronization obtained numerically with h = 10 -6 (see Eq. |?]) ), while the solid line is the stability boundary obtained 
through MSF. Parameters of the system are fi = 0.1; (i) ft = a = 0.1, (ii) /3 = 0.005, a = 0.114. 
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FIG. 4: Space-time- amplitude plot showing unstable synchronization in the ring of diffusively coupled oscillators, with fj, = a = 
/3 = 0.1. 
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(N,K)=(20,0.6) 




FIG. 5: Space-time- amplitude plot showing transitions from unstable to stable synchronization in a ring of diffusively coupled 
oscillators, with fj, = a = f3 = 0.1. 
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FIG. 6: Variation of the Kuramoto order parameter R versus the coupling strength K for a ring of diffusively coupled non- 
identical oscillators with parameters fi — a = /3 = 0.1, Atii = 0.05, n — 2. 




FIG. 7: Enlargement of the Kuramoto order parameter R versus the coupling strength K for the same parameters as in Fig. 
6 . 
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FIG. 8: Stability diagram of phase synchronization (R > 0.998) and disordered steady-state (R < 0.998J for a ring of diffusively 
coupled nonidentical oscillators. Parameters are fi = 0.1, Aiii„ = 0.05, k = 2, (i): a = f3 — 0.1, (ii): a — 0.114,/? = 0.005. 
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FIG. 9: Variation of the Master Stability Function (MSF) A versus K and N for a chain of open-ended of diffusive coupling 
oscillators, (i) fj, = f3 = a = 0.1 and (ii) /i = 0.1; /3 = 0.005, a — 0.114. 
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FIG. 10: Stability diagram in the (N, K) plane of chain of open ended diffusively coupled oscillators. The shaded area denotes 
the region of stable synchronization obtained numerically with h — 10 -6 (see Eq. ^ ), while the solid line is the stability 
boundary obtained through the Master Stability Function (MSF). Parameters of the system are n — 0.1; (i) (3 — a — 0.1, (ii) 
p = 0.005, a = 0.114. 




FIG. 11: Kuramoto order parameter R versus the coupling strength K for a chain of open-ended diffusively coupled non-identical 
oscillators with parameters jj, = a = j3 = 0.1, Aw = 0.05, k = 2. 
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FIG. 12: Stability diagram of phase synchronization (R > 0.998,) and disordered steady-state (R < 0.998,) for an open-ended 
chain of diffusively coupled nonidentical oscillators. Parameters are /i = 0.1, Aw — 0.05, k = 2, (i): a = f3 = 0.1, (ii): 
a = 0.114, fj = 0.005. 



